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High  order  essentially  non-oscillatory  (ENO)  schemes,  originally  designed  for  compress¬ 
ible  flow  and  in  general  for  hyperbolic  conservation  laws,  are  applied  to  incompressible 
Euler  and  Navier-Stokes  equations  with  periodic  boundary  conditions.  The  projection  to 
divergence-free  velocity  fields  is  achieved  by  fourth  order  central  differences  through  Fast 
Fourier  Transforms  (FFT)  and  a  mild  high-order  filtering.  The  objective  of  this  work  is  to 
assess  the  resolution  of  ENO  schemes  for  large  scale  features  of  the  flow  when  a  coarse  grid  is 
used  and  small  scale  features  of  the  flow,  such  as  shears  and  roll-ups,  are  not  fully  resolved. 
It  is  found  that  high-order  ENO  schemes  remain  stable  under  such  situations  and  quantities 
related  to  large-scale  features,  such  as  the  total  circulation  around  the  roll-up  region,  are 
adequately  resolved. 
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1  Introduction 


In  this  paper  we  consider  numerically  solving  the  incompressible  Navier-Stokes  or  Euler 
equations 


Ut  +  UUX  +  VUy  =  /i(UrX  +  Uyy)  ~ 

Vt+UVX  +  VVy  =  (l(VXX  +  Vyy)  ~  Py  (l.l) 

UX  +  Vy  =  0 

or  their  equivalent  conservative  form 

T  (n  )x  “I-  (tzu)y  —  p(^xx  T  ^yy )  Pi 

Vt  +  (uv)x  +  ( V2)y  =  p(vxx  +  Vyy)  ~  Py  (  1  .2) 

UX  +  Vy  --  0 

where  (u,  v)  is  the  velocity  vector,  p  is  the  pressure,  p  >  0  for  the  Navier-Stokes  equations  and 
p  =  0  for  the  Euler  equations.  The  numerical  methods  we  use  are  the  high-order  essentially 
non-oscillatory  (ENO)  schemes,  originally  designed  for  compressible  flow  and  in  general  for 
hyperbolic  conservation  laws  [8],  [12].  The  equation  is  defined  on  the  box  [0,27t]  x  [0,  27t] 
with  periodic  boundary  conditions  in  both  directions.  We  choose  two  space  dimensions  for 
easy  presentation,  although  our  method  is  also  applicable  for  three  space  dimensions. 

In  some  sense  equations  (1.1)  are  easier  to  solve  numerically  than  their  compressible 
counter-parts  because  the  latter  have  solutions  containing  possible  discontinuities  (for  ex¬ 
ample  shocks  and  contact  discontinuities).  However,  the  solution  to  (1.1),  even  if  for  most 
cases  smooth  mathematically,  may  evolve  rather  rapidly  with  time  t  and  may  easily  become 
too  complicated  to  be  fully  resolved  on  a  feasible  grid.  Traditional  linearly  stable  schemes, 
such  as  spectral  methods  and  high-order  central  difference  methods,  are  suitable  for  the 
cases  where  the  solution  can  be  fully  resolved,  but  typically  produce  signs  of  instability  such 
as  oscillations  when  small  scale  features  of  the  flow,  such  as  shears  and  roll-ups,  cannot  be 
adequately  resolved  on  the  computational  grid  (see,  for  example,  the  last  contour  in  Fig.  1). 
Although  in  principle  one  can  always  overcome  this  difficulty  by  refining  the  grid,  today’s 
computer  capacity  seriously  restricts  the  largest  possible  grid  size. 

In  the  last  few  years  there  is  a  lot  of  activity  in  designing  high-order,  nonlinearly  stable 
“shock  capturing”  schemes  for  compressible  flow  and  in  general  for  hyperbolic  conservation 
laws.  See,  for  example,  [2],  [7],  [8],  [12],  and  the  references  listed  therein.  The  philosophy  of 
such  schemes  is  to  give  up  fully  resolving  rapid  transition  regions  or  shocks,  just  to  “capture” 
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them  in  a  stable  and  somehow  globally  correct  fashion  (e.g.,  with  correct  shock  speed),  but  at 
the  same  time  to  require  a  high  resolution  for  the  smooth  part  of  the  flow.  The  success  of  such 
an  approach  for  the  compressible  flow  is  documented  by  many  examples  in  the  literature. 
One  example  is  the  one  and  two  dimensional  shock  interaction  with  vorticity  or  entropy 
waves  [12],  [13].  The  shock  is  captured  sharply  and  certain  key  quantities  related  to  the 
interaction  between  the  shock  and  the  smooth  part  of  the  flow,  such  as  the  amplification  and 
generation  factors  when  a  wave  passes  through  a  shock,  are  well  resolved.  Another  example  is 
the  homogeneous  turbulence  for  compressible  Navier-Stokes  equations  studied  in  [13].  In  one 
of  the  test  cases,  the  spectral  method  can  resolve  all  the  scales  using  a  2562  grid,  while  third 
order  ENO  with  just  642  points  can  adequately  resolve  certain  interesting  quantities  such  as 
the  average  fluctuation  Mach  number  and  maximum  divergence,  although  it  cannot  resolve 
local  quantities  achieved  inside  the  rapid  transition  region  such  as  the  minimum  divergence. 
The  conclusion  seems  to  be  that,  when  fully  resolving  the  flow  is  either  impossible  or  too 
costly,  a  “capturing”  scheme  such  as  ENO  can  be  used  on  a  coarse  grid  to  obtain  at  least 
some  partial  information  about  the  flow. 

In  this  paper  we  perform  a  similar  numerical  study  for  the  incompressible  equation  (1.1), 
to  see  what  one  can  get  by  using  the  high-order  ENO  schemes  on  a  coarse  grid,  without 
fully  resolving  the  flow.  We  choose  doubly  periodic  shear  layers,  used  in  [1],  as  our  test  case. 
A  spectral  method  with  5122  points  seems  able  to  fully  resolve  the  flow  up  to  t  —  8,  but 
begins  to  show  signs  of  underresolution  (wriggles  in  vorticity)  thereafter.  This  indicates  the 
tremendous  requirement  upon  computation  resources  if  one  tries  to  resolve  everything  in  the 
flow.  We  then  use  the  third  order  ENO  method  (which  is  fourth  order  in  the  L\  sense)  and 
coarse  grids  (642  and  1282  points),  to  assess  its  resolution.  We  find  that  the  scheme  remains 
stable  for  coarse  meshes  and  certain  quantities  related  to  the  smooth  part  of  the  flow,  such 
as  the  total  circulation  around  the  roll-up,  are  adequately  resolved  by  ENO  methods. 

A  pioneer  work  in  applying  shock  capturing  compressible  flow  techniques  to  incompress¬ 
ible  flow  is  by  Bell,  Colella  and  Glaz  [1],  in  which  they  considered  a  second  order  Godunov 
type  discretization,  investigated  the  projection  into  divergence-free  velocity  fields  for  general 
boundary  conditions,  and  discussed  accuracy  of  time  discretizations.  Since  our  objective  in 
this  paper  is  to  assess  the  resolution  of  the  ENO  method  for  (1.1),  we  choose  a  periodic 
boundary  condition  to  simplify  the  projection.  General  boundary  conditions  would  have  to 
be  handled  either  by  more  complicated  projection  [1]  or  by  artificial  compressibility  methods 
[4].  We  are  currently  investigating  ENO  schemes  for  such  cases. 
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2  The  ENO  Method 


We  solve  (1.2)  in  its  equivalent  projection  form 


where  P  is  the  Hodge  projection  into  divergence-free  fields,  i.e.,  if^*J=P^^J,  then 

ux  +  vy  =  0  and  vy  —  ux  =  vy  —  ux.  See,  e.g.,  [1].  For  the  current  periodic  case  the  additional 
condition  to  obtain  a  unique  projection  P  is  that  the  mean  values  of  u  and  v  are  preserved, 
fon  I?  «(*,  y)dxdy  =  /02t  u(x,  y)dxdy  and  /2*  f2*  v(x,  y)dxdy  =  f02w  /02’  v(x,  y)dxdy. 

We  use  Nx  and  Ny  (even  numbers)  equally  spaced  grid  points  in  x  and  y,  respectively. 
The  grid  sizes  are  denoted  by  Ax  =  ^  and  Ay  =  and  the  grid  points  are  denoted  by 
x,  =  iAx  and  y}  =  jAy.  The  approximated  numerical  values  of  u  and  v  at  the  grid  point 
(x,, yj)  are  denoted  by  utJ  and  . 

We  first  describe  the  numerical  implementation  of  the  projection  P.  In  the  periodic  case 
this  is  easily  achieved  in  the  Fourier  space.  We  first  expand  u  and  v  using  Fourier  collocation: 


uN(x,y)=  £  £ 


I{kx+ly) 


v/v(x,y)=  £  £ 


,I(kx+ly) 


/=_Al  k=-^ 


where  /  =  \f—\,  iiki  and  v^l  are  the  Fourier  collocation  coefficients  which  can  be  computed 
from  the  point  values  u;j  and  v,j,  using  either  FFT  or  matrix- vector  multiplications.  The  de¬ 
tail  can  be  found  in,  e.g.,  [3].  Derivatives,  either  by  spectral  method  or  by  central  differences, 
involve  only  multiplications  by  factors  dxk  or  df  in  (2.2)  because  e^kx+ly^  are  eigenfunctions 
of  such  derivative  operators.  For  example, 


for  spectral  derivatives; 


dxk  =  Ik,  df  =  II 


,  2/sin(i^)  M  2/sin(ifl) 

"  AJ  '  d>  -  Ay 


for  the  second  order  central  differences  which,  when  used  twice,  will  produce  the  second 
order  central  difference  approximation  w,il~Y+w>-'  for  Wxr,  and 


2/^(1  —  cos(A:Ax))(7  —  cos(A;Ax)) 
Ax 


2/ J(1  -  cos(/Ay))(7  -  cos(/Ay)) 


3 


for  the  fourth  order  central  differences  which,  when  used  twice,  will  produce  the  fourth  order 
central  difference  approximation  16(a''ti+W|-iM^2+u'‘-2)~30Wi  for  Wxr_  High  order  filters,  such 
as  the  exponential  filter  [10],  [9]: 


axk  =  e 


y  _  p-“(tB2p 


a,  =  e 


(2.6) 


where  2 p  is  the  order  of  the  filter  and  a  is  chosen  so  that  e~a  is  machine  zero,  can  be  used 
to  enhance  the  stability  while  keeping  at  least  2p-th  order  of  accuracy.  This  is  especially 
helpful  when  the  projection  P  is  used  for  the  under- resolved  coarse  grid  with  ENO  methods. 
We  use  the  fourth  order  projection  (2.5)  and  the  filter  (2.6)  with  2 p  —  8  in  our  calculations. 
This  will  guarantee  third  order  accuracy  (fourth  order  in  L\ )  of  the  ENO  scheme.  We  will 
denote  this  combination  (the  fourth  order  projection  plus  the  eighth  order  filtering)  by  P4. 

To  be  precise,  if  ^  ^  ^  =  P4  ^  ^ 

and  v ,  then  the  Fourier  collocation  coefficients  of  u  and  v  are  given  by 


and  Uki  and  v^i  are  Fourier  collocation  coefficients  of  u 


t  _  x  y<ff(<ffu-ffi)) 

u  k<T‘  (4)2  +  «)2 ' 


-  dxkv) 

kU‘  (dir  +  m2 


(2.7) 


where  and  erf  are  defined  by  (2.6)  with  2 p  =  8,  and  and  df  are  defined  by  (2.5). 

Next  we  shall  describe  the  ENO  scheme  for  (2.1).  Since  (2.1)  is  equivalent  to  the  non¬ 
conservative  form  (1.1),  it  is  natural  to  implement  upwinding  by  the  signs  of  u  and  t>,  and 
to  implement  ENO  equation  by  equation.  The  ?--th  order  ENO  approximation  of,  e.g.,  (u2)x 
is  thus  summarized  as  follows: 


1.  Take  f(x)  =  u2(x,y)  with  y  fixed.  We  start  with  the  point  values  /,  =  f{xi); 

2.  Define  a  function  h(x)  satisfying  f(x)  —  ^  fx_  aI  anc^  'ts  primitive  H(x)  = 

fx  h(^)d^.  For  each  j  +  a  (r  +  l)-th  order  polynomial  QJ+ 1/2(2:)  is  constructed 
interpolating  H(x)  at  points  near  £j+i/2.  It  is  remarked  in  [12]  that  the  Newton 
divided  differences  of  H(x)  are  constant  multiples  of  those  of  f(x)  of  one  order  lower. 
Since  we  do  not  need  the  zeroth  order  difference,  the  approximation  Q1+ 1/2  of  H(x) 
can  be  accomplished  using  the  information  of  /,  only,  without  explicitly  constructing 
H{x)  or  its  differences; 

3.  The  stencil  of  the  polynomial  Q(x)  is  determined  adaptively  by  upwinding  and  smooth¬ 
ness  of  f(x).  It  starts  with  either  Xj  or  xJ+\  according  to  whether  u  >  0  or  11  <  0,  then 
one  point  to  the  left  or  right  is  added  to  the  stencil  each  time  by  comparing  the  two 
relevant  divided  differences  and  picking  the  smaller  one  in  magnitude. 


4 


4.  The  derivative  f'(x)  is  finally  approximated  by  the  conservative  difference  ~(/J+ j/2  — 
fj-1/2)  where  the  numerical  flux  is  computed  by  fj+1/2  =  j^Q{x) |x=xJ+1/2- 

5.  In  the  actual  coding  of  the  algorithm,  undivided  differences  should  be  used  instead  of 
the  divided  differences  to  reduce  round-off  errors.  There  are  also  ways  to  make  the 
procedure  more  economical  on  a  vector  computer.  The  details  can  be  found  in  [13]. 


The  approximation  to  (uv)x  is  the  same  as  above  with  f(x)  =  a (x,y)v(x,y),  and  that 
for  {uv)y  and  for  (v2)y  are  similar,  with  upwinding  based  on  v. 

There  are  two  ways  to  handle  the  second  derivative  terms  for  the  Navier-Stokes  equations. 
One  can  absorb  them  into  the  convection  part  and  treat  them  using  ENO.  For  example, 
f(x)  =  u2(x,y)  can  be  replaced  by  f(x)  =  u2(x,y)  —  yu(x,y)x,  where  n(x,y)x  itself  can  be 
obtained  using  either  ENO  or  central  difference  of  a  suitable  order.  The  remaining  procedure 
for  computing  f{x)x  would  be  the  same  as  described  above.  Another  simpler  possibility  is 
just  to  use  standard  central  differences  (of  suitable  order)  to  compute  the  double  derivative 
terms.  Our  experience  with  compressible  flow  is  that  there  is  little  difference  between  the 
two  approaches,  especially  when  the  viscosity  y  is  small. 

In  the  above  we  have  described  the  discretization  for  the  spatial  derivatives 


y  =  yj 


Next  we  use  a  third  order  TVD  (total  variation  diminishing)  Runge-Kutta  method  [11]  to 
discretize  the  resulting  ODE: 


obtaining: 


P4  Lij 


(2.9) 


+  A  til 


fA< 


(2.10) 


Notice  that  we  have  used  the  property  P40  P4  =  in  obtaining  the  discretization  (2.10) 
from  (2.9). 
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This  explicit  time  discretization  is  expected  to  be  nonlinearly  stable  under  the  CFL 
condition 


At 


max 


+ 


+  2  fi 


+ 


1 


Ax2  Ay2 


<  1 


(2.11) 


Ax  Ay 

(see  [1 1]).  For  small  y  (which  is  the  case  we  are  interested  in)  this  is  not  a  serious  restriction 
on  At. 


3  Numerical  Results 

We  present  two  numerical  examples  in  this  section. 

Example  1:  This  example  is  used  to  check  the  third  order  accuracy  of  our  ENO  scheme 
for  smooth  solutions.  We  first  take  the  initial  condition  as 

u(x,  y,  0)  =  —  cos(x)  sin(y),  i>(x,  y,  0)  =  sin(x)  cos(y)  (3.1) 

which  was  used  in  [5].  The  exact  solution  for  this  case  is  known: 

u(x,  y,t)  —  —  cos(x)  sin (y)e~2)it,  v(x,  y,  t )  =  sin(x)  cos (y)e~2lit  (3-2) 

We  take  Ax  =  A y  —  jj  with  N  =  32,64,128  and  256.  The  solution  is  computed  up 
to  t  =  2  and  the  Li  error  and  numerical  order  of  accuracy  are  listed  in  Table  3.1.  For  the 
y  =  0.05  case,  we  list  results  both  with  fourth  order  central  approximation  to  the  double 

derivative  terms  (central)  and  with  ENO  to  handle  the  double  derivative  terms  by  absorbing 

them  into  the  convection  part  (ENO).  We  can  clearly  observe  fully  third  order  accuracy 
(actually  better  in  many  cases  because  the  spatial  ENO  is  fourth  order  in  the  L\  sense)  in 
this  table. 


Table  3.1:  Accuracy  of  ENO  Schemes  for  (3.1) 


N 

n  =  0 

y  =  0.05,  central 

y  =  0.05,  ENO 

Z/2  error 

order 

Li  error 

order 

Li  error 

order 

32 

9.10(-4) 

5.28(-4) 

4.87(-4) 

64 

5.73(-5) 

3.99 

3.20(-5) 

4.04 

3.09(-5) 

3.98 

128 

3.62(-6) 

3.98 

1.93(-6) 

4.05 

1.89(-6) 

4.03 

256 

2.28(-7) 

3.99 

1.18(-7) 

4.03 

1.16(-7) 

4.03 

We  then  take  the  initial  condition  as 
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u(x,t/,0)  =  —  sin2(x)  sin(2y),  v(x,y,  0)  =  sin(2x)  sin2(y)  (3.3; 

which  was  used  in  [1].  We  again  take  Ax  =  Ay  —  —  with  N  =  1 6, 32, 64, 128  and  256  and 
compute  the  solution  up  to  t  =  2.  However,  this  time  the  exact  solution  is  unknown.  As 
in  [1],  we  test  the  accuracy  by  computing  the  i2  difference  between  the  solutions  at  the 
grid  sizes  2Ax  and  Ax,  on  the  coarser  grid.  If  w/\x  =  w  +  CAxr  +  0(Axr+1),  then  this 
difference  w2 ax  —  w&T  =  (2r  —  l)CAxr  +  0(Axr+1)  would  predict  the  correct  asymptotic 
order  of  accuracy  r  and  the  error  itself  on  the  finer  grid  multiplied  by  2r  —  1.  The  result  is 
summarized  in  Table  3.2,  where  w  =  (u,e)r  and  L2  diff=  j|w2AX  —  Wax||l2-  In  this  table 
the  estimated  order  and  error  are  obtained  using  the  remarks  above.  For  p  —  0.05,  we  have 
again  provided  results  both  with  fourth  order  central  approximation  to  the  double  derivative 
terms  (central)  and  with  ENO  to  handle  the  double  derivative  terms  by  absorbing  them 
into  the  convection  part  (ENO).  The  result  in  Table  3.2  reconfirms  better  than  third  order 
accuracy  for  the  third  order  ENO  scheme. 

Table  3.2:  Accuracy  of  ENO  Schemes  for  (3.3) 


A 

(i  =  0 

p  =  0.05,  central 

p  =  0.05,  ENO 

L2  diff 

order 

error 

L2  diff 

order 

error 

L2  diff 

order 

error 

■ai 

1 . 14(- 1 ) 

3.20(-2) 

3.60(-2) 

■ai 

1.40(-2) 

3.02 

1.96(-3) 

Imrcittia 

3.52 

2.66(-4) 

2.93(-3) 

3.62 

2.60(-4) 

128 

1.46(-3) 

3.26 

1.69(-4) 

1 .8 1  (-4) 

3.94 

1.26(-5) 

1.80(-4) 

4.02 

1.18(-5) 

256 

1.11  (-4) 

3.77 

8.78(-6) 

1.09(-5) 

4.06 

6.91  (-7) 

1 . 1 0  ( -  5 ) 

4.04 

7 . 1 5  ( -  7 ) 

Example  2:  This  is  our  test  example  to  study  resolution  of  ENO  schemes  when  the  grid 
is  coarse.  It  is  a  double  shear  layer  taken  from  [1]: 


u(x,y,0)  =  | 


tanh((y  -  tt/2 )/p)  y  <  n 


v(x,  y,  0)  =  8  sin(x) 


(3.4) 


tanh((37r/2  -y)/p)  y  >  tt 
where  we  take  p  =  zr / 15  and  8  =  0.05.  The  Euler  equations  {p  =  0)  are  used  for  this  example. 
The  solution  quickly  develops  into  roll-ups  with  smaller  and  smaller  scales,  so  on  any  fixed 
grid  the  full  resolution  is  lost  eventually.  For  example,  the  expensive  run  we  performed  using 
5 1 22  points  for  the  spectral  collocation  code  (with  a  18-th  order  filter  (2.4))  seems  able  to 
resolve  the  solution  fully  up  to  t  =  8,  then  begins  to  lose  resolution  as  indicated  by  the 
wriggles  in  the  vorticity  contour  at  t  =  10  (the  last  contour  in  Fig.  1).  On  the  other  hand, 
the  ENO  runs  with  642  and  1282  points  produces  smooth,  stable  results  (Fig.  2  and  3). 
Apparently  with  these  coarse  grids  the  full  structure  of  the  roll-up  is  not  resolved.  However, 
when  we  compute  the  total  circulation 
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(3.o) 


cq=  u>(x,y)dxdy  =  /  udx  4-  vdy 
Jci  Jm 

around  the  roll-up  by  taking  Q  —  [j ,  rf]  x  [0,27 r]  and  using  the  rectangular  rule  (which  is 
infinite  order  accurate  for  the  periodic  case)  on  the  line  integrals  at  the  right-hand-side  of 
(.'3.5),  we  can  see  that  this  number  is  resolved  much  better  than  the  roll-up  itself  (see  Table 
3-3). 

Table  3.3:  Resolution  of  the  Total  Circulation 


t 

2 

4 

6 

8 

10 

ENO  642 

0.87300 

3.07100 

7.16889 

9.88063 

10.90122 

ENO  1282 

0.87452 

2.97810 

7.30999 

10.34414 

11.79418 

spectral  5 12 2 

0.87433 

2.98029 

7.28308 

10.46212 

11.85875 

4  Concluding  Remarks 

High  order  ENO  schemes  can  be  applied  to  the  incompressible  Euler  and  Navier-Stokes 
ecpiations  to  obtain  stable,  under-resolved  results  on  a  coarse  grid.  Some  global  quantities 
such  as  the  circulation  around  the  roll-up  region  can  be  faithfully  resolved. 
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